/************************************************************************
*************************************************************************
*************************************************************************

gendatva2.do

Code to generate a hypothetical dataset, with financial constraints and/or
variation in input prices.

12 August 2015: 		Created--borrows heavily from gendatva.do (@)
13 August 2015: 		Calibrated to Chile 311 parameters (@)
26 Octobe 2015: 		Updated to match gendat2
11 Novemb 2015:			Simplified choice function for M
						IMPORTANT fixed the nonlinear markov process


The production function is
Y = exp(om) X(K,L)

and unconstrained M is 
M = sqrt( 3*`normy' + `normy'^3  )

where is output normalized by its standard deviation 
normy = Y / SDev(Y)

(Clearly this choice has no theoretical basis...)

Arguments:
klfunc constraintmult rho pk pl pm N T sigeps sigom pricemsd

klfunc			-	a command that uses K and L to generate X(K,L) and stores it in a variable X_kl
constraintmult	-	the multiplier on "wealth" in the liquidity constraint
rho				-	persistence of productivity
N				-	number of firms
T				-	number of years in the dataset
pm pm2 			- 	The cobb douglas coefficient, or if mfunc=quad the coefficient on the squared term
sigeps			-	standard deviation of the random shock
sigom			-	standard deviation of the persistent shock
pricemsd		-	standard deviation of the price of flexible inputs



*************************************************************************
*************************************************************************
************************************************************************/


cap program drop gendatva2
program gendatva2

	

	//User will pass the multiplier on the constraint, the parameters, and the sample size
	args klfunc constraintmult rho N T sigeps sigeta pricemsd nlprod fracobs
	
	
	
	if mi("`fracobs'") {
		local fracobs = 1/2
	}
	
	
	//////////////////////
	//////Set up firms and productivity
	
	//Variance of productivity
	local sigom = sqrt( 1 / (1-`rho'^2) ) * `sigeta'
	
	//Initial state
	clear
	set obs `N'
	gen id 			= _n
	gen omega 		= rnormal()*`sigom'
	gen k 			= 2.115286*rnormal() 	+ 6.946116 	
	gen l 			= 0.9518106*rnormal() 	+ 3.659235	
	gen logpricem 	= 0
	
	//Add time dimension
	expand `=`T'+2'
	bysort id: gen t = _n
	
	xtset id t
	
	//////////////////////
	//////Create the flow of productivity and factors
	
	//Random shock and its mean
	if !mi("`sigeps'") {
		gen eps 		= rnormal()*`sigeps'
	}
	else {
		gen eps 		= rnormal()
	}
	gen tmp 		= exp(eps)
	quietly: sum tmp
	gen meaneps 	= r(mean)
	drop tmp
	
	//Generate input price shocks
	if !mi("`pricemsd'") {
		gen priceeps 		= rnormal()*`pricemsd'
	}
	else {
		gen priceeps 		= rnormal()*0
	}
	
	//Generate future capital, labor, and productivity
	//gen rainshock 		= rnormal()*`sigeta'											if t > 1
	//gen eta 			= .5*rainshock + rnormal()*`sigeta'
	
	gen eta 				= rnormal()*`sigeta'
	
	if mi("`nlprod'") {
		//replace omega 		= `rho'*l.omega + (1-`rho')*eta 							if t > 1
		replace omega 		= `rho'*l.omega + eta 										if t > 1
	}
	else {
		replace omega 		= `nlprod' + eta											if t > 1
	}
	
	/*
	These constants are calibrated from Industry 311 in the Chilean ENIA
	*/
	
	local kk1 = 1.003592
	local ke1 = .0359078
	local ke2 = .1957049
	local kc  = -.012955
	
	local ll1 = .9519013
	local le1 = .4676349
	local le2 = .1845552
	local lc  = .188374
	
	
	replace k 			= `kc' + `kk1'*l.k + `ke1'*l.eta  + `ke2'*l2.eta 	if t > 2
	replace l 			= `lc' + `ll1'*l.l + `le1'*l.eta  + `le2'*l2.eta	if t > 2
	
	
	//Persistence of prices (calibrated from data)
	local persistp .9274344
	
	local pricescale = sqrt(1-`persistp'^2)											//Rescale error to keep the variance stationary
	replace logpricem	= `persistp'*l.logpricem + `pricescale'*priceeps	if t > 1
	
	gen K 				= exp(k)
	gen L 				= exp(l)
	
	
	
	/*
	These constants are calibrated from Industry 311 in the Chilean ENIA
	(I treat log[stassets] as S )
	*/
	
	local savr1			= .4920887
	local savr2 		= .3205638
	local savc 			= 1.940491
	local savsd 		= 2.438074 			//Total standard deviation
	
	local savres 		= 1.354485 			//Standard deviation of the residual

	
	gen savings 		= 7.417812 + rnormal()* `savsd'
	replace savings 	= `savc' + `savr1'*l.savings + `savr2'*l2.savings + rnormal()* `savres' if t > 2
	gen unobs 			= 7.417812 + rnormal()* `savsd'
	replace unobs 		= `savc' + `savr1'*l.unobs + `savr2'*l2.unobs + rnormal()* `savres' if t > 2
	gen wealth 			= exp(`fracobs'*savings + (1-`fracobs')*unobs )
	
	//Rescale by means
	quietly: sum K
	local scK = r(mean)
	
	quietly: sum wealth
	local scW = r(mean)
	
	
	gen collateral		= (K/`scK')^.5 *  (wealth/`scW')^.5
	gen con				= `constraintmult' * collateral
	
	//Price of intermediates
	gen pricem			= exp(logpricem)
	
	
	//Contribution of capital and labor
	`klfunc'
	
	//////////////////////
	//////Choose flexible input and compute output
	
	/*
	tempvar normy   
	gen `normy' 		= exp(omega) * X_kl 
	sum `normy'
	replace `normy' = `normy' / r(sd)
	gen M 			= min( 3*sqrt( 3*`normy' + `normy'^3  )  , con / pricem )
	*/
	
	
	/*
	
	//Parameters of choice function
	local bm 		= -27.58009
	local km		= -.1902844
	local lm 		= 1.350338
	local llm		= -.0427895
	local kkm		= .0505756
	local klm		= -.0472306
	local om		= 4.480906
	
	gen mopt		= `bm' + `km'*k + `lm'*l + `kkm'*k^2 + `klm'*k*l + `llm'*l^2 + `om'*omega
	*/
	
	
	local bm 		= 0 //-27.58009
	local km		= -1.00901
	local lm 		= 5.527479
	local llm		= -.0433079
	local kkm		= .0468967
	local klm		= -.039854
	local om		= -21.16055
	local oom		= 1.918072
	local okm		= .1173157
	local olm		= -.5962474
	
	gen mopt		= `bm' + `km'*k + `lm'*l + `kkm'*k^2 + `klm'*k*l + `llm'*l^2 + `om'*omega + `oom'*omega^2 + `okm'*omega*k + `olm'*omega*l

	
	
	gen M 			= min( exp( mopt ) , con / pricem )
	gen Y 			= exp(omega) * X_kl * exp(eps)
	
	
	foreach var of varlist M Y {
		local tmp = lower("`var'")
		gen `tmp' = log(`var')
	}
	
	//Throw away the initial periods
	drop if t <= 2
	
	
	//Is the firm constrained?
	gen constrained = abs(con / pricem - M) < .0001
	order y k l m, last
	
	
	

end



